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ABSTRACT 

j  Status  and  some  recent  developments  in  the  techniques  for  assessing  the  reliability 
of  finite  element  solutions  are  summarized.  Discussion  focuses  on  a  number  of  aspects 
including:  the  major  types  of  errors  in  the  finite  element  solutions;  techniques  used  for 
a  posteriori  error  estimation  and  the  reliability  of  these  estimators;  the  feedback  and 
adaptive  strategies  for  improving  the  finite  element  solutions;  and  postprocessing 
approaches  used  for  improving  the  accuracy  of  stresses  and  other  important  engineering 
data.  Also,  future  directions  for  research  needed  to  make  error  estimation  and  adaptive 
improvement  practical  are  identified.  (~ 

1.  INTRODUCTION 

The  finite  element  method  has  become  the  main  tool  in  computational  mechanics 
and  to  date  the  method  is  used  for  solving  a  large  class  of  engineering  problems  which  are 
stated  in  terms  of  differential,  pseudo-differential,  integral  or  integro-differential 
equations.  The  success  of  the  finite  element  method  is  manifested  by  the  development  of 
over  five  hundred  user-oriented  finite  element  program  systems,  and  over  two  hundred 
pre-  and  postprocessing  packages  (Ref.l).  There  are  over  20,000  finite  element  users 
worldwide  who  are  estimated  to  spend  about  $500  million  annually  on  finite  element 
analyses.  The  literature  on  the  subject  is  nearly  overwhelming  and  to  date  there  are  over 
two  hundred  monographs  and  conference  proceedings  published  on  various  aspects  of  the 
finite  element  technology  (Ref.  2).  A  review  of  some  recent  developments  in  finite 
element  method  is  contained  in  survey  papers  (see,  for  example,  Refs.  3  and  4)  and  a 
state-of-the-art  monograph  (Ref.  5).  Despite  the  significant  advances  made  on  the  theory 
and  algorithmic  tools  of  the  finite  element  method,  the  selection  of  the  finite  element 


model  for  a  particular  problem  is  largely  based  on  intuition  and  experience  gained  from 
solving  similar  problems.  Moreover,  the  assessment  of  the  reliability  of  the  finite  element 
solution  continues  to  be  the  most  difficult  aspect  of  the  finite  element  analysis.  Some 
aspects  of  the  errors  in  finite  element  computations  are  discussed  in  Refs.  6  and  7.  Error 
estimation  and  control  is  now  at  the  "cutting  edge"  of  the  finite  element  technology. 

In  recent  years  considerable  interest  has  been  shown  in  the  development  of  reliable 
error  estimates  as  well  as  feedback  procedures  (or  adaptive  strategies)  by  which  a  required 
accuracy  of  the  finite  element  solution  can  be  most  economically  reached.  The  number  of 
publications  on  the  subject  has  been  steadily  increasing  and,  although  two  symposia 
proceedings  have  been  published  on  the  subject  (Refs.  8  and  9),  there  is  a  need  to  broaden 
awareness  among  practicing  engineers  and  research  workers  about  the  status  and  recent 
developments  in  various  aspects  of  quality  assessment  and  control  of  the  finite  element 
solutions.  This  paper  is  an  attempt  to  fill  this  void.  Specifically,  the  objectives  of  this 
paper  are:  a)  to  review  and  assess  the  current  techniques  used  for  error  estimation  and 
adaptive  improvement  of  the  finite  element  solutions;  and  b)  to  identify  future  directions 
of  research  needed  to  make  error  estimation  and  adaptive  improvement  practical. 

The  contents  of  the  paper  are  arranged  as  follows:  in  section  2,  the  major  types  of 
errors  in  a  finite  element  solution  are  listed;  in  section  3,  the  available  technology  for 
estimating  the  discretization  errors  is  reviewed.  Then  the  different  strategies  for 
adaptive  improvement  of  the  finite  element  solution  are  discussed  in  section  4.  In  section 
5  the  postprocessing  techniques  used  for  improving  the  quality  of  stresses  and  other 
important  engineering  data  are  reviewed.  The  future  directions  of  research  on  quality 
assessment  and  control  are  identified  in  section  6. 

There  is  a  large  literature  available.  The  cited  references  are  selected  for 
illustrating  the  ideas  presented  and  are  not  necessarily  the  only  significant  contributions 
to  the  subject.  The  discussion  in  this  paper  is  kept  on  a  descriptive  level  and  for  all  the 
mathematical  details  the  reader  is  referred  to  the  cited  literature. 


2.  MAJOR  ERRORS  IN  THE  FINITE  ELEMENT  SOLUTION 

In  general,  the  reliability  of  the  finite  element  solutions  to  engineering  problems  is 
influenced  by  the  following  factors  (see  Refs.  6  and  7): 

a)  Reliability  of  the  mathematical  model  which  describes  the  real  structure  in 
mathematical  terms.  The  stochastic  formulation  is  one  of  the  possible  mathematical 
formulations  of  the  engineering  problem.  However,  deterministic  formulations  are  more 
commonly  used. 

b)  Errors  and  uncertainties  in  the  input  information  (of  the  mathematical  model). 
These  include  uncertainties  in  material,  geometry;  boundary  conditions;  loading 
environment;  and  the  characteristics  of  the  probability  fields  in  the  stochastic 
formulation. 

c)  Discretization  errors  which  are  caused  by  the  numerical  discretization  of  the 
continuous  mathematical  model,  as  well  as  by  the  truncation  of  infinite  processes  (e.g., 
iterative  procedures  and  summation  of  infinite  series). 

d)  Roundoff  errors  which  occur  in  the  implementation  of  numerical  algorithms  on 
computers  with  finite  precision  (limitation  in  representing  real  numbers  due  to  the 
finiteness  of  the  computer  word  length). 

An  assessment  of  the  reliability  of  the  mathematical  model  requires  the 
identification  of  the  range  of  validity  of  the  mathematical  theory  used  in  describing  the 
model,  and  the  effect  of  violating  some  of  the  basic  assumptions  of  the  theory  on  the 
response  predictions.  Examples  of  these  issues  are  provided  by  the  range  of  validity  of 
dimensional  reduction  when  applied  to  plates  with  rapidly  varying  thickness  (i.e.,  the 
reduction  from  three-dimensional  continuum  theory  to  two-dimensional  plate  theory  - 
Refs.  10  and  11);  the  applicability  of  the  linear  elasticity  theory  in  the  neighborhood  of  a 
crack  tip;  and  the  effect  of  "pseudo-supports"  on  the  response  of  structural  members  (Ref. 


The  stochastic  formulation,  which  is  widely  used  in  dynamics,  leads  to  a 
deterministic  formulation  for  the  mean  and  correlation  functions  of  the  solution  (see  Refs. 
12  to  20).  Monte  Carlo  methods  are  sometimes  valuable  tools  for  the  mathematical 
treatment  of  the  uncertainties  in  the  problem  parameters.  As  an  alternative  to  stochastic 
modeling,  an  assessment  of  the  effect  of  the  uncertainties  in  the  input  information  on  the 
response  can  be  made  by  evaluating  the  sensitivity  of  the  response  of  the  system  to  each 
of  the  input  parameters. 

The  discretization  error  represents  the  difference  between  the  exact  and  the 
approximate  (finite  element)  solutions  for  the  mathematical  model.  The  determination  of 
the  actual  discretization  error  requires  the  availability  of  the  exact  solution,  which  is 
rarely  known  in  practical  problems.  However,  even  when  the  exact  solution  is  not 
available,  it  is  possible  to  construct  quantitative  estimates  for  the  discretization  error  and 
to  determine  the  rate  of  change  of  the  error  as  the  number  of  degrees  of  freedom  in  the 
finite  element  model  increases. 

Since  discretization  errors  are  functions,  error  norms  are  usually  used  to  measure 
the  size  of  these  errors.  For  a  brief  discussion  of  the  error  norms  used  in  finite  element 
analyses,  see  Appendix  I.  The  class  of  discretization  errors  also  includes  the  errors  caused 
by  the  truncation  of  iterative  procedures  and  the  summation  of  infinite  series.  To  date 
discretization  error  analyses  are  of  a  priori  and  a  posteriori  types.  This  will  be  discussed 
further  in  the  succeeding  section. 

In  recent  years  considerable  work  has  been  devoted  to  the  control  of  round-off  errors 
and  its  accurate  estimation.  An  overview  of  the  current  state-of-the-art  is  contained  in 
Ref.  21,  in  which  practical  methods  and  algorithms  are  presented  for  carrying  out  basic 
processes  of  numerical  analysis  with  guaranteed  error  bounds.  The  error  bounds  account 
for  rounding  errors,  and  wherever  appropriate,  discretization  errors.  Some  of  the 
algorithms  presented  in  Ref.  21  have  been  incorporated  in  the  high-accuracy  arithmetic 
package  (ACRITH)  marketed  by  IBM  (Ref.  22).  An  important  feature  of  the  package  is  its 


ability  to  accumulate  inner  products  exactly.  Some  flaws  in  the  package  have  been 
discussed  in  Ref.  23.  Analytical  approaches  for  developing  a  priori  and  a  posteriori  bounds 
for  roundoff  errors  are  presented  in  the  two  monographs  (Refs.  24  and  25),  and  in  the  two 
papers  (Refs.  26  and  27).  A  priori  and  a  posteriori  estimates  of  roundoff  errors,  based  on 
linearized  perturbation  theory,  are  given  in  Ref.  28.  The  primary  use  of  these  estimates  is 
to  assess  the  suitability  of  numerical  algorithms  based  on  the  sensitivity  of  the  solution  to 
perturbations  in  the  input  data.  An  extensive  bibliography  on  error  analysis  based  on  the 
use  of  interval  mathematics  is  contained  in  Ref.  29. 

3.  ESTIMATION  OF  DISCRETIZATION  ERRORS 

Since  the  development  of  modern  finite  element  method  in  the  mid  1950's,  attempts 
have  been  made  by  engineers  and  researchers  to  obtain  information  about  discretization 
errors,  and  to  prove  the  convergence  of  the  solutions.  Mathematicians  became  interested 
in  the  subject  in  mid  1960's.  Until  recently,  the  two  groups  worked  almost  independently, 
and  even  now,  the  emphasis  of  the  work  of  the  two  groups  is  different.  The  availability  of 
information  about  discretization  errors,  in  addition  to  its  importance  in  assessing  the 
reliability  of  the  finite  element  solutions,  provides  the  possibility  of  obtaining  the  "best 
possible"  solution  within  an  allowable  cost  range. 

3.1  Types  of  Discretization  Error  Estimates 

In  general,  there  are  two  types  of  discretization  error  estimates:  a)  a  priori 
estimates;  and  b)  a  posteriori  estimates. 

A  priori  estimates  are  based  on  a  knowledge  of  the  characteristics  of  the  solution 
such  as  its  smoothness,  and  provide  qualitative  information  about  the  asymptotic  rate  of 
convergence  as  the  number  of  degrees  of  freedom  goes  to  infinity.  Although  a  priori 
estimates  are  accurate  for  the  worst  case  in  a  particular  class  of  solutions  of  the  problem 
(e.g.,  class  of  all  solutions  which  have  the  second  derivative  bounded  by  a  constant),  they 
usually  do  not  provide  information  about  the  actual  error  for  a  given  finite  element 
model.  Nevertheless,  a  priori  estimates  can  be  very  effective  when  used  in  extrapolation 


techniques.  To  date  a  priori  error  estimates  (or  convergence  rate  estimates)  are  available 
for  a  wide  variety  of  finite  elements.  For  a  brief  discussion  of  a  priori  error  estimates  see 
Appendix  I.  A  thorough  analysis  of  these  errors  is  given  in  many  monographs  and 
publications  (see,  for  example,  Refs.  30  to  35). 

A  posteriori  estimates  use  information  obtained  during  the  solution  process,  in 
addition  to  some  a  priori  assumptions  about  the  solution.  For  an  abstract  treatment  of  the 
questions  related  to  incomplete  information  see  Ref.  36.  Discussion  in  the  succeeding 
subsections  will  focus  on  a  posteriori  estimates  which  can  provide  quantitatively  accurate 
measures  of  the  discretization  error. 


3.2  Discretization  Error  Diagnostics 

Early  work  by  engineers  on  discretization  errors  focused  on  measuring  the  degree  of 
reliability  of  the  solution  obtained  with  a  given  finite  element  model,  and  developing  error 
diagnostics  through  checking  the  satisfaction  of:  a)  the  natural  boundary  conditions  (for 
variational-based  finite  elements);  b)  the  equilibrium  equations  by  calculating  the 
unbalanced  forces  (equilibrium  defects)  both  in  the  interior  of  the  element  as  well  as  on  its 
interfaces,  see  Ref.  37  (for  compatible  displacement  models);  and  c)  the  compatibility 
equations  by  calculating  the  strain  incompatibility  in  the  interior  of  the  elements  as  well 
as  across  adjacent  elements  (for  equilibrium,  mixed  and  nonconforming  models). 

The  aforementioned  error  diagnostics  are  based  on  considering  the  finite  element 
solution  as  an  exact  solution  of  the  original  problem  with  perturbed  input  data  (e.g., 
magnitude  of  unbalanced  forces).  The  main  problem  is  to  relate  the  magnitude  of  the 
perturbations  to  the  response  quantities  of  interest  (e.g.,  stresses).  Modern  a  posteriori 
error  estimation  techniques  are  often  indirectly  related  to  error  diagnostics  through  the 
effects  of  superconvergence.  Superconvergence,  which  can  occur  at  special  points,  refers 
to  extraordinarily  high  rates  of  convergence  (and  concomitant  high  accuracy)  observed  at 
these  points. 


Equilibrium  defects  may  be  deceptive  unless  properly  used  and  interpreted  (see  Ref. 
38).  The  proper  interpretation  of  equilibrium  defects  is  still,  in  general,  an  open 
question.  For  a  discussion  of  some  aspects  of  this  subject  see  Refs.  37,  39  and  40. 

3.3  Classical  Approaches  for  Estimating  Discretization  Errors 

The  two  classical  approaches  for  estimating  the  discretization  errors  are: 
a)  Extension  Methods  based  on  comparing  the  finite  element  solutions  obtained  using 
a  sequence  of  meshes  of  increasing  refinement  (h  extension);  with  an  increasing  degree  of 
the  polynomial  shape  functions  (p  extension);  or  using  a  combination  of  the  two  (h-p 
extension).  The  computer  implementation  of  the  h-,  p-  and  h-p  extensions  are  referred  to 
as  the  h-,  p-  and  h-p  versions,  respectively.  Extension  methods  are  used  to  determine  the 
rate  of  convergence  of  the  finite  element  solution,  and  its  actual  error. 

The  h-version  of  the  finite  element  method  is  the  standard  one.  The  mathematical 
theories  for  the  p-  and  h-p  versions  are  given  in  Refs.  41  to  51.  Engineering  and 
computational  aspects  of  extension  methods  are  discussed  in  Refs.  52  to  56.  The 
a  posteriori  error  bounds  in  Refs.  39,  46,  52,  54  and  56  are  obtained  by  using  extrapolation 
lrocedures.  The  selection  of  the  correct  extrapolation  procedure  is  based  on  the 
theoretical  a  priori  error  estimates,  and  is  essential  for  the  effectiveness  of  the 
a  posteriori  estimates  (see  Refs.  39,  46,  52,  54  and  56). 

The  effectiveness  of  extension  methods  in  practical  applications  depends  to  a  great 
extent  on  their  efficient  computer  implementation.  Herein  efficient  implementation 
refers  to  the  combination  of  the  coarse  and  fine  grid  solutions  in  the  h-  version  or  the 
lower-  and  higher-degree  polynomial  solutions  in  the  p-version  to  generate  the  final 
solutions.  For  the  h-version  efficient  computer  implementation  can  be  achieved  by  using 
multigrid  concepts  which  are  widely  used  in  conjunction  with  the  finite  difference  method 
(see  Ref.  57).  For  the  p-  and  h-p  versions,  efficient  implementation  can  be  made  by  using 
elements  with  nested  set  of  shape  functions  (i.e.,  hierarchic  shape  functions,  see  Ref. 

52).  A  brief  discussion  of  the  hierarchic  shape  functions  and  their  computational 


advantages  is  given  in  Appendix  II. 


b)  Dual  (or  Complementary)  Procedure  -  based  on  obtaining  two  solutions  using 
direct  and  complementary  variational  methods  to  provide  bounds  on  the  error  of  the  global 
response  characteristics.  This  approach  can  be  applied  to  self-adjoint  problems 
formulated  using  minimum  (rather  than  stationary)  variational  principles.  For  elasticity 
problems  this  is  accomplished  by  using  the  direct  and  dual  formulations  of  the  problem 
(based  on  the  principles  of  the  minimum  potential  energy  and  the  minimum  complementary 
energy)  to  obtain  upper  and  lower  bounds  of  the  error  measured  in  the  energy  norm  (see 
Ref.  58).  Applications  of  the  dual  procedure  are  given  in  Refs.  59,  60  and  61.  The 
procedure  is  expensive  but  leads  to  guaranteed  upper  bounds  of  the  error  (measured  in  the 
energy  norm)  which  are  not  pessimistic. 

3.4  A  Posteriori  Error  Estimates  and  Error  Indicators 

In  recent  years  considerable  effort  has  been  devoted  to  the  development  of 
computable  (i.e.,  element  by  element)  a  posteriori  error  estimates  based  on  local  error 
estimators.  These  error  estimators  are  based  on  information  obtained  during  the  solution 
process  itself.  The  work  presented  in  Refs.  62  to  66  started  this  activity  and  was  extended 
in  various  directions  (se«  Refs.  67  to  77).  For  example,  in  Ref.  67  error  estimates  are 
presented  for  linear  elements  (with  linear  shape  functions)  using  a  priori  bounds  on  single 
elements  as  well  as  a  posteriori  information  about  the  second  derivatives.  All  the  error 
estimates  developed  in  the  cited  references  are  for  compatible  (or  conforming)  finite 
element  models. 

The  finite  element  solutions  have  equilibrium  defects  in  the  interior  and  on  the 
portion  of  the  boundary  where  tractions  are  prescribed;  as  well  as  jumps  in  the  tractions  at 
interelement  boundaries.  The  equilibrium  defects  are  nearly  self-equilibrating.  The 
square  of  the  energy  norm  of  the  total  error  (i.e.,  the  strain  energy  of  the  error)  can  be 
obtained  by  summing  the  squares  of  the  energy  norms  of  the  errors  for  individual 
elements,  each  obtained  by  solving  (or  estimating  the  solution  of)  an  auxiliary  problem 


characterizing  the  equilibrium  defects  in  that  element  (see  Refs.  62,  69,  73  and  77).  In 
case  of  odd-degree  elements  (elements  with  odd-degree  polynomial  shape  functions),  the 
boundary  defects  dominate  the  residual,  and  for  even-degree  elements  the  interior  defects 
are  dominant  (see  Ref.  78).  The  energy  norms  of  the  aforementioned  auxiliary  problem 
can  be  accurately  estimated  or  obtained  by  a  finite  element  with  properly  selected  shape 
functions.  Hierarchic  shape  functions  are  particularly  useful  for  this  purpose  (see  Refs. 

68,  71,  79  and  80). 

The  estimation  of  the  square  of  the  energy  norm  of  the  total  error  through 
summation  of  elemental  contributions  can  also  be  applied  to  nonlinear  problems  (including 
limit  and  bifurcation  points  in  the  solution  space),  see  Refs.  62,  64,  79,  81  and  82. 

The  local  error  estimators  (computed  for  individual  elements)  are  referred  to  in  Ref. 
66  as  error  indicators  n  (  Cl  ),  where  Q.  is  the  element  domain.  Some  approaches 
for  constructing  error  indicators  are  described  in  Refs.  79  and  83.  Error  indicators  can  be 
used  to  construct  overall  problem  error  estimators,  e  ,  for  the  problem.  The  expression 
of  c  in  terms  of  n  depends  on  the  particular  error  norm  being  used  (see  Refs.  66,  67,  71, 
73,  76,  83,  84  and  85).  For  example,  the  energy  norm  error  is  given 

^  2  V  2  (e) 

r  =  2j  n  (£2  ) 

elements 

Other  error  estimators  proposed  include  residual  errors  in  the  constitutive  relations 
(Refs.  40,  74  and  86);  and  in  the  strain-displacement  relations  (Ref.  87).  Error  estimators 
for  singular  perturbation  problems  are  given  in  Refs.  83  and  88. 

3.5  Reliability  of  Error  Estimators 

The  a  posteriori  error  estimators  described  in  the  preceding  subsections  are 
asymptotically  correct  if  the  ratio  of  the  estimator  to  the  true  error  converges  to  one  as 
the  true  error  tends  to  zero  (when  the  mesh  size  tends  to  zero,  or  when  the  degree  of 
polynomial  shape  function  goes  to  infinity).  Therefore,  in  order  to  assess  the  reliability  of 
these  estimators  for  a  finite  grid,  an  effectivity  index,  9  ,  was  introduced  in  Refs.  65  and 


89.  The  effectivity  index  is  defined  as  follows: 

9  =7Pri 

where  e  is  the  a  posteriori  error  estimator  and  ||  e  ||  is  the  norm  of  the  error.  The 
estimator  is  asymptotically  correct  when  9-*-l  as  II  e  ||  -*-  0. 

For  practical  applications  it  is  essential  that  1 6  -1  f  be  small  (e.g.,  less  than  0.2 
when  the  error  |]  e  ||  is  of  the  order  of  10%  of  the  norm  of  the  solution,  and  should 
decrease  as  the  error  decreases,  see  Ref.  62). 

It  is  preferable  that  0  be  greater  than  one  so  that  the  true  error  is  overestimated. 
The  asymptotic  correctness  of  the  error  estimator  is  related  to  the  superconvergence 
effects.  In  addition,  some  assumptions  about  the  meshes  and  the  solution  are  needed, 
without  which  the  quality  of  the  estimators  can  be  low.  The  assumptions  about  the  meshes 
are  satisfied  when  they  are  constructed  adaptively.  In  such  cases  I  6  -l|  is  of  the  order  of 
0.1.  The  mathematical  theory  of  error  estimators  is  discussed  in  Refs.  73,  77,  78  and  85. 
Application  to  one-dimensional  problems  are  given  in  Refs.  50  and  63. 

3.6  Error  Estimators  for  Stationary  and  Transient  Problems 

Most  of  the  error  estimates  reported  in  the  literature  are  for  stationary  problems 
(linear  and  nonlinear).  Error  estimates  have  also  been  developed  for  some  transient 
problems  (see  Refs.  79,  83  and  90  to  97).  Stationary  problems  include  static  stress  analysis 
and  are  mostly  governed  by  elliptic  systems  of  differential  equations.  The  transient 
problems  include  dynamic  response  and  wave  propagation  problems  and  are  governed  by 
equations  of  the  parabolic  and  hyperbolic  types. 

Among  the  elliptic  structural  mechanics  problems  for  which  a  posteriori  error 
estimates  exist  are  one-dimensional  rod  with  variable  stiffness,  second-order  partial 
differential  equations,  and  two-dimensional  elasticity.  Numerical  results  have  been 
presented  for  plane  stress  problem  of  a  cracked  plate,  plane  strain  problem  of  a  dam,  L- 
shaped  region  with  a  corner  singularity,  see  Refs.  62,  68-70,  84-86  and  88.  The  cited 
references  contain  numerical  examples  that  demonstrate  the  effectiveness  of  the  error 
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estimators.  Most  of  the  problems  presented  in  the  cited  references  are  governed  by 
second-order  differential  equations  and  are  modeled  by  using  conforming  elements.  Error 
estimates  for  plate  bending  problems  and  for  the  boundary  element  method  are  discussed 
in  Ref.  98.  Only  few  finite  element  programs  have  facilities  for 

error  estimation,  based  on  the  computation  of  error  indicators  (see,  for  example,  Refs.  81 
99  and  100).  This  will  be  discussed  further  in  Section  4.3. 

3.7  Comments  on  Error  Estimators 

1.  The  selection  of  an  error  measure  depends  on  the  goal  of  the  computation,  and  on 
which  response  quantities  need  to  be  accurately  determined,  e.g.,  buckling  loads,  stresses 
in  a  critical  zone;  stress  intensity  factors  at  crack  tips;  or  density  in  compressible  flow 
problems  (for  a  discussion  of  the  goals  of  stress  analysis  see  Ref.  101).  This  is  particularly 
important  since  the  smallness  of  one  measure  does  not  guarantee  the  smallness  of  the 
errors  in  other  response  quantities,  computed  from  the  finite  element  solution.  An 
example  of  this  was  given  in  Ref.  39  where  a  small  energy  norm  error  was  associated  with 
large  errors  in  the  computed  stresses. 

2.  For  the  elliptic  problems  considered  in  the  literature  the  computation  of  the 
error  indicators  takes  a  small  percentage  of  the  total  solution  time.  Significant  dif¬ 
ferences  can  occur  in  the  cost  of  computing  different  error  estimators.  Also,  error 
estimators  for  higher-order  elements  are  more  expensive  to  compute  than  for  lower -order 
elements. 

3.  Pointwise  estimation  of  the  errors  in  detailed  response  quantities  (e.g.,  stresses, 
displacements  and  stress  intensity  factors)  is  considerably  more  difficult  than  the  evalu¬ 
ation  of  the  local  energy  norm  errors.  For  simple  problems  (e.g.,  Saint-Venant  torsion 
problem,  one-dimensional  problems,  and  two-dimensional  plane  elasticity  problems) , 
attempts  have  been  made  to  obtain  pointwise  error  estimates  for  stresses.  These  error 
estimates  are  based  on  a  postprocessing  approach  (Ref.  102)  or  on  heuristic  arguments 
using  the  experience  gained  from  finite  difference  discretizations  (Refs.  72  and  103). 


4.  When  the  goal  of  the  analysis  is  the  accurate  determination  of  detailed  response 
quantities  such  as  the  stresses,  pointwise  error  estimates  are  needed.  If  these  pointwise 
estimates  are  not  available,  then  energy  norm  errors  are  computed,  and  additional 
numerical  tests  are  performed  such  as  element-by-element  equilibrium  test  and  action- 
reaction  test  (Ref.  104). 

4.  THE  ADAPTIVE  IMPROVEMENT  OF  FINITE  ELEMENT  SOLUTIONS 

Adaptive  improvement  of  finite  element  solutions  refers  to  improving  the  quality  of 
the  solutions  by  enriching  the  approximation  in  some  manner  so  as  to  achieve  the  "best" 
solution  for  a  given  computational  effort  (or  cost).  Adaptive  procedures  have  in  common 
with  feedback  approaches  the  fact  that  the  computed  information  is  used  for  steering  the 
process.  However,  it  is  the  attempt  to  obtain  the  "best"  solution  that  distinguishes 
adaptive  processes  from  feedback  approaches.  Therefore,  an  adaptive  process  is  a 
feedback  approach  which  is  optimal  with  respect  to  clearly  defined  objectives  (see  Refs. 
89,  105  and  106).  The  adaptive  process  is  normally  performed  after  an  initial  solution  is 
already  available;  and  regions  of  the  solution  domain  where  the  accuracy  is  not 
satisfactory  have  been  identified  according  to  a  preselected  set  of  criteria.  This  process 
can  be  iterated  using  the  last  solution  as  the  new  initial  solution.  The  term  adaptive 
refers  to  the  fact  that  each  model  modification  step  uses,  in  some  optimal  way,  the 
information  provided  by  previous  steps.  For  a  discussion  of  the  notion  of  adaptivity  in 
numerical  procedures  see  Refs.  89  and  105. 

Adaptive  processes  can  result  in  considerable  savings  in  computational  effort  if  an 
initial  coarse  grid  is  used  followed  by  an  effective  strategy  to  improve  the  quality  of  the 
solution  wherever  needed.  In  most  of  the  literature  on  feedback  approaches,  neither  the 
goals  of  the  process  nor  the  criteria  used  to  test  whether  the  feedback  approach  is 
adaptive,  are  defined.  A  feedback  approach  can  be  adaptive  (and  hence  optimal)  with 
respect  to  a  set  of  criteria  and  nonadaptive  with  respect  to  ether  criteria.  In  the  h- 
version  of  the  finite  element  method,  the  feedback  consists  of  the  construction  of 


successive  meshes  to  achieve  the  goal  of  computation,  and  in  the  h-p  version  both  the 
shape  functions  and  the  mesh  are  enriched  simultaneously. 


Mesh  selection  and  improvement  can  be  carried  out  by  using  Al-based  expert 
systems,  or  a  combination  of  expert  systems  and  a  feedback  approach.  The  expert  systems 
act  as  an  expert  consultant  which  decides  the  direction  to  steer  the  adaptive  process  (Ref. 
107),  or  can  be  used  as  a  preprocessor  for  designing  the  finite  element  model. 

In  general,  a  feedback  process  has  a  number  of  key  elements  including  (see  Refs.  76, 
84,  88,  89  and  105): 

a)  error  measure.  Each  model  modification  step  should  aim  to  minimize  some  error 
measure.  This  error  measure  is  related  to  the  goal  of  the  computation; 

b)  the  feedback  procedure  (or  adaptive  modification  strategy)  which  aims  at 
minimizing  the  error  measure  selected  in  a)  in  the  most  effective  manner;  and, 

c)  stopping  criterion.  This  criterion  is  related  to  the  goal  of  the  computation, 
usually  through  reducing  the  error  measure  to  a  prescribed  tolerance. 

The  efficiency  of  the  feedback  process  is  measured  through  a  well-defined  cost 
function  (see  Ref.  106).  If  either  the  error  measure  or  the  cost  function  is  changed,  the 
feedback  algorithm  and  the  resulting  finite  element  model  may  change  significantly.  As 
an  example  to  this,  the  use  of  the  error  measures  based  on  the  energy  norm  and  the 
maximum  stress  norm  can  lead  to  two  quite  different  finite  element  models. 

4.1  Strategies  for  Adaptive  Improvement  of  Solutions  (Feedback  Procedures) 

Among  the  different  strategies  used  for  adaptive  improvement  of  the  finite  element 
solutions  (feedback  procedures)  are  the  following  four: 

a)  successive  selection  of  the  meshes; 

b)  moving  the  nodes  (node  relocation); 

c)  successive  selection  of  the  local  order  of  the  approximation;  and, 

d)  simultaneous  selection  of  the  meshes  and  the  local  order  of  approximation. 
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In  the  first  strategy,  the  degree  of  the  elements  is  kept  constant  and  the  mesh  is 
constructed  in  a  feedback  manner.  It  is  usually  referred  to  as  the  h-method.  It  can  result 
in  refining  or  de-refining  the  mesh.  De-refinement  of  the  mesh  may  be  necessary  in 
transient  (time-dependent)  problems  and  in  nonlinear  problems  when  continuation  method 
is  used  for  their  solution.  A  theoretical  analysis  of  the  various  feedback  approaches  and 
their  optimality  in  one-dimensional  problems  is  given  in  Refs.  80  and  108.  The  feedback 
approaches  for  multidimensional  problems  are  based  on  heuristic  arguments  and  analogies 
with  one-dimensional  problems.  Note  that  the  feedback  algorithms  involve  a  direct  or 
indirect  use  of  the  error  indicators  discussed  in  Section  3.4;  and  the  adaptive  strategies 
usually  aim  at  making  the  error  indicators  for  different  elements  nearly  equal. 

In  the  second  strategy,  the  quality  of  the  finite  element  solution  is  improved  by 
optimizing  the  location  of  the  nodes  keeping  the  number  of  degrees  of  freedom  fixed  and 
the  degree  of  the  elements  constant.  This  is  usually  referred  to  as  R-method  (Refs.  109 
and  110).  For  hyperbolic  equations  the  R-method  is  more  often  used  in  conjunction  with 
finite-difference  methods  than  with  finite  element  methods.  Because  of  the  kinship 
between  the  two  methods,  the  basic  ideas  of  adaptive  approaches  used  in  conjunction  with 
finite  differences  can  be  applied  to  the  finite  element  methods.  For  a  discussion  of 
adaptive  finite -difference  methods  see,  for  example,  Refs.  Ill  to  116. 

In  the  application  of  the  moving  element  method  to  initial  value  (transient) 
problems,  a  fixed  number  of  finite  elements  are  moved  so  as  to  concentrate  the 
computational  grid  in  regions  containing  nonuniformities  in  the  solution  (e.g.,  shocks,  near 
shocks,  boundary  layers,  and  sharp  moving  fronts).  Both  the  nodal  amplitudes  and  nodal 
positions  move  continuously  with  time  in  such  a  way  as  to  satisfy  simultaneous  ordinary 
differential  equations  (in  time)  which  minimize  the  partial  differential  equation 
residuals.  Galerkin-type  finite  elements  are  typically  used  for  the  spatial  discretization 
and  stiff  ordinary  differential  equation  solvers  are  used  for  the  temporal  integration. 
Several  mesh  moving  techniques  have  been  developed  for  one-dimensional  parabolic  and 
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hyperbolic  systems  (see,  for  example,  Refs.  95,  97,  and  117  to  122).  For  two-dimensional 
systems,  substantial  difficulties  still  remain  including  mesh  tangling  and  gross  distortion 
(see,  for  example,  Refs.  79,  94  and  123  to  125). 

The  third  strategy  is  based  on  successive  selection  of  the  order  of  the  polynomial 
approximation  inside  the  elements.  It  is  usually  referred  to  as  the  p-method  and  has  the 
advantages  over  the  first  two  methods  of  being  easy  to  implement  and  of  providing  a 
simple  formula  for  the  error  indicator.  This  is  particularly  true  when  hierarchic  basis 
functions  are  used,  since  the  stiffness  matrices  and  load  vectors  corresponding  to  a 
polynomial  of  degree  p  are  embedded  in  the  stiffness  matrices  and  load  vectors 
corresponding  to  polynomials  of  degrees  >  p+1  (see  Refs.  79,  83  and  126  to  128). 
Application  of  this  strategy  to  the  solution  of  three-dimensional  thermoelastic  problems  of 
anisotropic  solids  is  given  in  Ref.  129. 

The  fourth  strategy,  based  on  simultaneous  selection  of  the  meshes  and  the  local 
order  of  approximation  is  referred  to  the  h-p  method.  The  mathematical  foundations  for 
the  adaptive  h-p  method  and  R-methods  are  much  less  developed  than  for  the  other 
methods. 

A  variety  of  feedback  approaches  are  now  available  for  differential  equations  of 
different  types:  elliptic,  parabolic,  and  hyperbolic;  linear  and  nonlinear;  describing  a 
broad  spectrum  of  engineering  problems.  Significant  contributions  were  made  in  the  last 
few  years  (see,  for  example,  Refs.  67  to  71,  77,  79,  83,  84,  90  to  92,  94  to  97,  116  to  121, 
123,  125  and  130  to  132). 

Some  of  the  feedback  approaches  are  based  on  older  ideas  such  as  the  local 
extrapolation  techniques  and  the  deferred  correction  approach.  The  deferred  correction 
approach  is  used  in  various  software  packages  for  solution  of  initial  and  two-point 
boundary  value  problems.  It  is  based  on  using  the  numerical  solution  obtained  to  construct 
a  pseudo  or  neighboring  problem  whose  exact  solution  is  known  (e.g.,  polynomial  or  spline 
interpolation  of  the  discrete  numerical  solution).  The  pseudo  problem  is  then  solved  using 
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the  same  finite  element  model  as  that  used  for  the  original  problem.  The  error  in  the 
pseudo  problem  is  assumed  to  be  a  close  approximation  of  the  error  in  the  original  problem 
and  is  used  as  a  correction  to  that  solution.  The  technique  has  been  successfully  applied  to 
the  numerical  solution  of  stiff  systems  of  ordinary  differential  equations  (Refs.  133  and 
134)  and  has  recently  been  extended  to  finite  element  boundary-value  problems  (Refs.  135 
to  137). 

4.2  Optimal  Finite  Element  Mesh 

An  optimal  finite  element  mesh  is  the  one,  among  the  class  of  admissible  meshes,  for 
which  the  distribution  of  the  nodes  (and  associated  degrees  of  freedom)  minimize  the  error 
(measured  in  a  particular  norm).  In  the  class  of  admissible  meshes,  among  other  things, 
the  local  order  of  approximation  is  specified.  Obviously,  changing  the  error  norm  or  the 
class  of  admissible  meshes  may  change  significantly  the  optimal  mesh. 

The  first  attempts  to  produce  optimal  meshes  for  variational-based  finite  elements 
were  based  on  considering  the  nodal  positions  to  be  unknowns,  and  minimizing  the 
functional  (potential  energy)  with  respect  to  both  the  nodal  degrees  of  freedom  and  the 
nodal  positions  (Refs.  138  to  141).  Such  an  approach  was  found  to  be  impractical  since  it 
results  in  nonlinear  equations,  even  for  linear  problems. 

The  construction  of  an  optimal  mesh,  which  minimizes  the  error,  results  in  an  ill- 
posed  problem.  A  small  movement  of  a  node  of  the  optimal  mesh  leads  to  a  second-order 
change  in  the  error.  To  alleviate  this  problem  nearly  optimal  meshes  are  considered.  For 
transient  problems  the  feedback  procedure  used  for  generating  nearly  optimal  meshes,  is 
based  on  using  various  spring  constants  and  other  stabilization  parameters  in  the  finite 
element  model  (see  Refs.  94  and  117). 

More  recently,  it  was  proposed  to  characterize  an  optimal  mesh  as  one  for  which  the 
error  indicators  for  different  elements  are  nearly  equal.  For  one-dimensional  problems  it 
was  shown  that  the  mesh  with  equal  error  indicators  is  asymptotically  optimal  in  the  sense 
that  it  leads  to  the  minimum  error  for  a  given  number  of  degrees  of  freedom  (see  Refs.  51 
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and  65).  Most  of  the  currently-used  adaptive  strategies  attempt  to  construct  nearly 
optimal  meshes  (with  nearly  equal  error  indicators).  This  is  referred  to  as  the 
equilibration  principle  for  mesh  construction  (Ref.  66).  Important  applications  of  this 
principle  for  the  cases  when  the  solution  has  singularities  are  given  in  Refs.  53,  142  and 
143.  A  detailed  analysis  of  the  optimal  meshes  and  distribution  of  element  degrees  for  a 
one-dimensional  problem  is  given  in  Ref.  51.  In  the  cited  reference  the  characteristics  of 
the  optimal  meshes,  as  well  as  the  upper  and  lower  bounds  of  the  error  are  given,  for  both 
the  h-  and  h-p  methods. 

The  mathematical  foundations  for  the  construction  of  optimal  meshes  is  much  better 
developed  for  elliptic  problems  than  for  parabolic  and  hyperbolic  problems.  For  initial 
value  (transient)  problems  the  reader  is  referred  to  Refs.  79,  83,  90  to  92,  125  and  131.  In 
principle,  if  the  mesh  constructed  in  the  feedback  mode  has  nearly  the  same  accuracy 
(measured  in  some  norm)  as  that  of  the  optimal  mesh,  then  the  feedback  process  is  defined 
to  be  optimal.  Hence,  the  lower  estimates  for  the  error  of  an  optimal  mesh  have  close 
relation  to  the  construction  of  an  adaptive  approach. 

4.3  Computer  Implementation  of  Feedback  Procedures  (Adaptive  Strategies) 

One  of  the  most  difficult  aspects  of  feedback  procedures  is  the  computer 
implementation  of  the  automated  enrichment  or  adaptive  algorithm.  This  includes  the 
selection  of  the  data  structure  used  for  the  representation  (and  refinement)  of  the  grid, 
and  grid  management.  The  implementations  of  feedback  processes  are  different  for 
elliptic,  parabolic  and  hyperbolic  problems  (see  Refs.  75,  81,  93,  96,  99,  100,  144  to  147). 
Most  of  the  aforementioned  papers  on  adaptive  improvement  of  numerical  solutions  give 
numerical  examples  to  demonstrate  the  proposed  strategies,  but  do  not  present  general 
computer  programs  for  implementing  these  strategies. 

To  the  authors'  knowledge  none  of  the  large-scale,  general-purpose  commercial 
finite  element  systems  has  facilities  for  adaptive  enrichment  of  solutions.  However,  some 
special-purpose  and  pilot  finite  element  systems  have  implemented  adaptive  algorithms 
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for  practical  performance  studies.  Among  these  systems  are  the  Finite  Element  Adaptive 
Research  Solver  (FEARS)  developed  at  the  University  of  Maryland  in  1979-1980  (Ref.  99); 
the  EXPDES  system  developed  in  Belgium  (Ref.  100);  the  PLTMG  system  at  the  University 
of  California,  San  Diego  (Ref.  81);  the  self-adaptive  finite  element  simulator  (SAFES  code) 
developed  at  the  University  of  Wyoming  (Ref.  148);  and  the  PROBE  system  developed  by 
NOETIC  Technologies  Corporation  in  St.  Louis  (Ref.  52).  The  first  two  systems  and  the 
last  system  are  designed  to  solve  two-dimensional  elliptic  problems.  The  three  systems 
(FEARS,  PLTMG  and  EXPDES)  are  based  on  the  h-version.  The  last  system  (PROBE)  is  a 
commercial  program  which  uses  the  p-  and  h-p  extensions.  Although  it  is  currently  not 
self-adaptive,  the  developer  is  adding  the  automatic  enrichment  capability  to  it.  The 
SAFES  code  performs  the  spatial  refinement  for  elliptic  and  parabolic  linear  systems.  All 
of  the  aforementioned  programs  were  developed  for  sequential  computers.  It  is 
anticipated  that  modern  parallel  architecture  will  have  a  strong  impact  on  the 
implementation  of  adaptive  methods  (see  Refs.  145  and  149). 

4.4  Comments  on  Adaptive  Strategies  and  Feedback  Procedures 

The  following  comments  can  be  made  concerning  the  adaptive  strategies  reported  in 
the  literature: 

a)  Adaptive  processes  have  been  used  extensively  in  the  solution  of  ordinary 
differential  equations  (initial  value  and  two-point  boundary  value  problems).  Modern 
program  packages  incorporate  feedback  approaches. 

b)  Detailed  rigorous  mathematical  analyses  of  adaptive  algorithms  are  only 
available  for  some  one-dimensional  problems.  For  multidimensional  problems  the 
algorithms  used  for  the  adaptive  construction  of  the  desired  sequence  of  nearly 
equilibrated  meshes  (with  nearly  uniform  error  indicators)  are  largely  of  a  heuristic 
nature. 

c)  The  efficient  computer  implementation  of  adaptive  processes  has  not  matured 
yet.  However,  based  on  the  experience  acquired  so  far  certain  directions  are  evolving. 
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For  example,  in  the  adaptive  h-version  the  use  of  a  tree  structure  for  the  data  was  found 
to  provide  the  most  efficient  data  management  system.  In  the  adaptive  p-  and  h-p 
versions  the  use  of  hierarchic  shape  functions  is  optimal  because  the  stiffness  and  load 
matrices  for  a  certain  level  of  refinement  remain  unchanged  when  a  higher  level  is 
introduced,  and  iterative  methods  can  be  efficiently  used  for  solving  the  algebraic 
equations  associated  with  different  levels  of  refinement  (similar  to  the  multigrid  finite 
difference  setting);  see,  for  example,  Refs.  57,  75,  81,  96,  99,  100,  144  and  147. 

d)  For  elliptic  problems  adaptive  strategies  based  on  mesh  enrichment  (h-,  p-  and  h- 
p  methods)  are  more  commonly  used  than  those  based  on  node  relocation  (R-method).  For 
a  large  class  of  engineering  problems,  the  h-p  extension  has  an  exponential  rate  of 
convergence  with  respect  to  the  number  of  degrees  of  freedom,  while  the  h-  and  p- 
extensions  have  only  a  polynomial  rate  of  convergence  (see  Appendix  I).  However,  if 
properly-graded  meshes  are  used,  the  performance  of  the  p-extension  in  the  pre- 
asymptotic  range  is  very  close  to  the  best  performance  attainable  with  the  h-p  extension. 
Moreover,  for  quasi-uniform  meshes  the  rate  of  convergence  of  the  p-extension  cannot  be 
worse  than  that  of  the  corresponding  h-extension  (based  on  the  same  total  number  of 
degrees  of  freedom).  In  the  presence  of  singularities,  the  rate  of  convergence  of  the  p- 
extension  can  be  twice  that  of  the  h-extension  (Refs.  44  and  46  to  48). 

e)  For  initial  value  (parabolic  and  hyperbolic)  problems  moving  meshes  have  to  be 
used  either  continually  or  by  adaptively  determined  intervals  (Ref.  92).  The  major 
drawback  of  moving  mesh  techniques  is  that  they  do  not  add  and/or  delete  elements  as  the 
computation  progresses,  and  are  generally  not  capable  of  generating  a  solution  with  a 
prescribed  level  of  accuracy.  To  remedy  this  drawback  a  strategy  combining  mesh  moving 
and  refinement  (or  de-refinement)  has  recently  been  proposed  (see  Refs.  115,  124  and  146). 

f)  The  p-version  results  in  more  dense  matrices  than  the  h-version  and  is  more 
convenient  and  simpler  to  implement.  In  addition,  the  treatment  of  curved  elements 
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require  special  care  when  higher-degree  polynomials  are  used  (viz.,  smooth  mapping  such 
as  that  provided  by  the  blending  function  method  -  see  Appendix  II). 

g)  The  p-version  was  found  to  be  not  too  sensitive  to  the  "locking"  problems 
encountered  in  thin  plate  and  curved  shell  structures  (see  Ref.  150). 

5.  POSTPROCESSING  OF  FINITE  ELEMENT  SOLUTIONS  TO  IMPROVE  THEIR  QUALITY 

Often  the  primary  objective  of  the  finite  element  analysis  is  to  determine  certain 
response  quantities,  such  as  stress  intensity  factor,  stresses,  strains,  displacements,  or  flux 
in  some  area  of  the  domain  (or  structure).  The  simplest  approach  is  to  compute  this  data 
directly  from  the  finite  element  solution.  For  example,  stresses  are  calculated  by 
differentiating  the  displacement  field  and  using  the  constitutive  relations.  This  results  in 
lower  accuracy  for  the  stresses  (and  strains)  than  for  the  displacements.  Recent  studies 
have  demonstrated  that  the  accuracy  and  rate  of  convergence  of  stresses  (and  strains) 
depend  on  how  (and  where)  they  are  computed.  The  number  of  publications  on 
postprocessing  approaches,  based  on  the  superconvegence  concepts  (i.e.,  increased 
accuracy  and  improved  rates  of  convergence),  is  steadily  increasing  (for  example,  Ref.  151 
lists  200  publications  on  superconvergence).  Most  of  the  publications  deal  with  elliptic 
problems.  Among  the  proposed  postprocessing  approaches  for  elliptic  problems  are: 

a)  evaluating  the  stresses  at  numerical  quadrature  points  and  determining  their 
values  at  the  displacement  nodes  by  extrapolation; 

b)  averaging  or  smoothing  based  on  projection  techniques; 

c)  using  the  extraction  techniques; 

d)  computing  the  stresses  using  the  discarded  structural  equations  (corresponding  to 
prescribed  displacement  boundary  conditions);  and, 

e)  using  iterative  refinement  techniques  which  build  continuous  stress  field  and 
displacement  field  from  the  classical  solution. 

The  first  approach  is  by  far  the  most  commonly  used.  Superconvergence  has  been 
observed  for  displacement  derivatives  and  stresses  evaluated  at  the  quadrature  points  of 


uniform  meshes  (using  Gaussian  quadrature  formula  of  the  proper  degree-  see  Refs.  152  to 
159). 

Averaging  and  smoothing  techniques  can  be  very  effective  for  postprocessing  of  the 
finite  element  solution.  Several  approaches  for  averaging  and  smoothing  have  been 
proposed  in  the  literature  (see,  for  example,  Refs.  160  and  161).  The  analysis  of  simple 
averaging  techniques  is  discussed  in  Refs.  151  and  162.  The  superconvergence  for  general 
(nonuniform)  meshes  seems  to  be  an  open  question  (see  Ref.  163). 

Extraction  techniques  are  based  on  using  analytic  expressions  of  the  function  which 
approximate  the  kernel  of  the  functional  of  interest.  Because  of  its  higher  cost  of 
implementation,  this  approach  is  feasible  for  computing  the  stress  intensity  factor,  or  for 
obtaining  accurate  stresses  only  in  a  small  area  (critical  zone)  of  the  structure  (see  Refs. 
52,  102  and  164  to  169).  Extraction  methods  for  computing  the  stress  intensity  factor  in 
problems  with  singularities  at  the  corner  and  along  the  edges  are  presented  in  Ref.  170. 

The  fourth  approach  has  recently  been  extended  to  the  computation  of  stresses  at 
the  interfaces  between  elements  using  the  previously  computed  stiffness  matrices  and  load 
vectors  (Ref.  171).  For  one-dimensional  problems,  the  technique  requires  little 
computation  and  its  implementation  is  straightforward.  However,  for  two-dimensional 
problems,  the  stress  calculation  involves  the  evaluation  of  contour  integrals  and  the 
implementation  is  more  complex.  Moreover,  new  difficulties  arise  from  the  presence  of 
corners  in  polygonal  domains. 

Iterative  refinement  techniques,  such  as  those  reported  in  Refs.  172  to  174,  have 
been  quite  effective  in  reducing  equilibrium  errors  and  significantly  improving  the 
accuracy  of  the  calculated  stresses.  These  techniques  were  recently  shown  to  be 
equivalent  to  mixed  formulations  in  which  the  stresses  (or  strains)  and  displacements  are 
used  as  primary  variables  (Ref.  175).  The  techniques  have  high  potential  for  nonlinear 
problems  since  the  iterative  improvement  can  be  done  simultaneously  with  those 
necessitated  by  the  nonlinearity,  and  no  additional  computational  cost  is  involved. 
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The  aforementioned  approaches  were  discussed  in  the  context  of  elliptic  problems. 
Less  work  has  been  done  on  superconvergence  concepts  and  postprocessing  approaches  for 
parabolic  problems,  and  much  less  for  hyperbolic  problems  (see  Refs.  176  and  177). 

6.  FUTURE  DIRECTIONS  FOR  RESEARCH 

Considerable  attention  has  been  devoted  to  the  estimation  and  control  of 
discretization  errors  which  is  manifested  by  the  development  of  a  mathematical  theory  for 
a  posteriori  error  estimates  and  feedback  approaches  (adaptive  refinement  strategies). 
However,  the  theory  is  very  incomplete  and  has  been  numerically  tested  only  on  a  simple 
set  of  problems.  Moreover,  adaptive  strategies  have  only  been  incorporated  into  special- 
purpose  and  pilot  finite  element  systems  used  for  practical  performance  studies.  To-date 
none  of  the  general-purpose  commercial  finite  element  systems  have  facilities  for 
estimating  the  error  in  the  finite  element  solution,  or  for  adaptive  improvement  of  this 
solution.  To  remedy  this  situation  major  advances  are  needed  in  the  theory,  strategies  and 
algorithms  for  implementation  of  error  estimation  and  control.  Some  of  these  advances 
are  listed  subsequently. 

1.  Development  of  reliable  measures  for  estimating  the  errors  in  the  finite  element 
predictions  of  the  major  response  quantities  (selected  by  the  analyst  according  to  the  goal 
of  computation).  These  errors  are  due  to  the  simplifying  assumptions  made  in  abstracting 
the  mathematical  model  from  the  real  system  (structure);  uncertainties  in  the  input 
information  (of  the  mathematical  model);  and  numerical  discretization  of  the  continuous 
mathematical  model.  The  following  observations  can  be  made  about  estimating  the 
different  types  of  errors: 

a)  The  use  of  hierarchic  modeling  strategies  and  formulations  (e.g.,  hierarchy  of 
two-dimensional  plate  and  shell  theories)  allows  the  estimation  of  modeling  errors  as  well 
as  the  adaptive  refinement  of  the  model  whenever  needed. 

b)  The  use  of  stochastic  formulations  allows  the  study  of  the  effect  of  uncertainties 
in  input  information  on  the  response  predictions  of  the  finite  element  model. 


c)  Practical  a  posteriori  error  estimators  should  satisfy  the  following  four  criteria. 

i)  provide  reliable  local  assessment  of  the  error; 

ii)  are  computationally  inexpensive  to  evaluate; 

iii)  applicable  to  a  wide  class  of  discrete  finite  elements;  and, 

iv)  easy  to  use  in  conjunction  with  adaptive  improvement. 

Note  that  the  choice  of  the  error  estimator  depends  on  the  goal  of  the  computation  (e.g., 
buckling  analysis,  detailed  stress  analysis,  determination  of  stress  intensity  factor,  etc.). 

In  the  absence  of  practical  error  estimators,  attempt  should  be  made  to  develop 
computationally  inexpensive  error  sensors  or  correction  indicators  which  give  an  indication 
of  the  regions  where  mesh  refinement  is  needed.  The  sensitivity  derivatives  of  the 
solution  to  the  degree  of  approximation  can  be  used  as  error  sensors.  Also,  numerical 
checks  to  validate  the  solution  should  be  conducted.  Examples  of  these  checks  are 
element-by-element  equilibrium  check  and  action-reaction  test. 

2.  Development  of  efficient  adaptive  improvement  strategies.  In  particular,  the 
development  of  finite  element  models  based  on  the  p-hierarchic  shape  functions  for 
higher-order  problems  (e.g.,  shear-flexible  plates  and  shells)  along  with  efficient  numerical 
algorithms  for  reanalysis.  The  element  computation  should  include  the  esstential 
information  needed  for  a  posteriori  error  estimation. 

3.  Efficient  computer  implementation  of  adaptive  strategies.  This  includes  using 
novel  computer  science  concepts  for  data  management. 

4.  Development  of  adaptive  strategies  for  new  computing  systems.  A  number  of 
attempts  have  been  made  to  develop  parallel  numerical  algorithms  and  parallel  software 
for  adaptive  processes  (see,  for  example,  Refs.  145  and  149).  However,  the  effective  use 
of  parallelism  in  adaptive  strategies  requires  the  use  of: 

a)  primitive-variable  formulations  (e.g.,  three  field  mixed  formulation  in  structural 
mechanics); 

b)  domain  decomposition  or  substructuring  (with  minimization  of  interfaces);  and 
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c)  operator  splitting  techniques  to  uncouple  the  resulting  algebraic  equations. 

5.  A  systematic  assessment  of  the  postprocessing  techniques  used  for  improving  the 
accuracy  of  derivative  calculation  (stresses  and  strains  in  displacement  finite  element 
models)  is  needed.  Also,  the  possibility  of  using  these  techniques  for  error  sensing  should 
be  investigated. 

6.  Development  of  Al-knowledge  based  system  for  efficient  use  of  finite  element 
programs  by  the  nonexperts. 

7.  Selection  of  meaningful  benchmark  problems  which  have  the  essential  features  of 
practical  problems  to  test  the  theory  and  the  effectiveness  of  the  adaptive  strategies 
developed. 

7.  CONCLUDING  REMARKS 

Status  and  some  recent  developments  in  the  techniques  for  error  estimation  and 
adaptive  improvement  of  finite  element  solutions  are  summarized.  Discussion  focuses  on 
a  number  of  aspects  including:  the  major  types  of  errors  in  finite  element  solutions; 
techniques  used  for  a  posteriori  error  estimation  and  the  reliability  of  these  estimators; 
the  feedback  and  adaptive  strategies  used  for  improving  the  finite  element  solutions;  and 
postprocessing  approaches  used  for  improving  the  accuracy  of  stresses  and  other  important 
engineering  data. 

The  status  of  finite  element  modeling,  error  estimation  and  adaptive  improvement  of 
finite  element  solutions  can  be  summarized  in  the  following: 

1.  Mathematical  modeling  errors.  Despite  their  importance,  have  not  received 
enough  attention  in  the  literature.  To  date  no  simple  and  general  approach  is  available  for 
assessing  these  errors. 

2.  Selection  of  the  formulation  and  initial  finite  element  models.  These  are  largely 
based  on  intuition  and  experience  gained  from  similar  problems.  The  modeling  process  can 
be  aided  by  the  information  available  on  the  asymptotic  rates  of  convergence  (a  priori 
estimates  of  the  discretization  errors)  for  different  finite  element  models  based  on 
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different  formulations.  These  a  priori  error  estimates  predict  the  asymptotic  rate  of 
convergence  as  the  mesh  size  tends  to  zero  (or  the  degree  of  polynomial  shape  function 
goes  to  infinity).  However,  except  in  the  adaptive  improvement  of  finite  element 
solutions,  a  priori  error  estimates  are  currently  not  used  in  constructing  the  initial  finite 
element  models.  The  knowledge-based  systems,  in  spite  of  their  high  potential  have  not 
been  sufficiently  developed  to  aid  in  the  modeling  process. 

3.  Assessment  of  the  reliability  of  the  finite  element  solution.  This  includes: 

a)  Selection  of  an  appropriate  measure  for  the  discretization  error,  which  depends  on 
the  goal  of  the  computation,  and  on  which  response  quantities  need  to  be  accurately 
determined;  and, 

b)  Calculation  of  a  posteriori  estimates  for  the  discretization  error  using  locally 
computable  (element-by-element)  error  indicators. 

The  most  commonly  used  error  measures  are  the  interior  and  boundary  residuals;  and 
the  local  energy  norm  error.  The  first  represents  the  equilibrium  defects  in  the  interior; 
on  the  portion  of  the  boundary  where  tractions  are  prescribed,  as  well  as  the  jumps  in  the 
tractions  at  interelement  boundaries.  The  energy  norm  error  is  defined  as  the  square  root 
of  the  strain  energy  of  the  error.  Pointwise  error  estimates  for  detailed  response 
quantities  (e.g.,  stresses  and  displacements)  are  available  only  for  simple  problems  (e.g., 
Saint  Venant  torsion  problem  and  one-dimensional  problems). 

Most  of  the  a  posteriori  error  estimates  reported  in  the  literature  are  for  stationary, 
(elliptic)  problems  (e.g.,  plane  stress  analysis  of  structures  subjected  to  conservative 
loading  and  modeled  by  conforming  displacement  finite  elements).  No  a  posteriori  error 
estimates  are  available  for  nonconforming  elements. 

4.  Postprocessing  of  finite  element  solutions.  The  accuracy  of  strains,  stresses  and 
stress  intensity  factors  computed  from  compatible  finite  element  models  can  be  improved 
by  postprocessing.  Several  approaches  have  been  proposed  in  the  literature.  A  systematic 
and  detailed  assessment  of  these  approaches  is  needed. 


5.  Quality  control  for  finite  element  solution.  This  is  accomplished  by  using 
adaptive  strategies  which  are  feedback  processes  aimed  at  achieving  the  desired  accuracy 
with  the  least  computational  effort  and/or  cost.  The  adaptive  process  is  performed  after 
an  initial  solution  is  already  available  and  the  regions  of  the  solution  domain  where 
accuracy  is  not  satisfactory  have  been  identified  (according  to  the  preselected  error 
measure).  Four  adaptive  strategies  have  been  proposed  in  the  literature.  These  include: 
i)  successive  selection  of  the  meshes  (h-method);  ii)  moving  the  nodes  (R-method);  iii) 
successive  selection  of  the  order  of  the  polynomial  approximation  inside  some  elements  (p- 
method);  and  iv)  simultaneous  selection  of  the  meshes  and  of  the  local  order  of 
approximation  (h-p  method). 

For  stationary  (elliptic)  problems  the  h-,  p-  and  h-p  methods  are  more  commonly 
used  than  the  R-method.  For  time-dependent  (parabolic  and  hyperbolic)  problems  adaptive 
strategies  are  based  on  either  the  moving  mesh  techniques  or  combination  of  mesh  moving 
and  refinement  (or  de-refinement). 

The  mathematical  theory  for  a  posteriori  error  estimates  and  feedback  approaches  is 
very  incomplete  and  has  been  numerically  tested  only  on  a  simple  set  of  problems.  In  spite 
of  the  considerable  attention  devoted  by  engineers  and  mathematicians  to  the  subject  of 
error  estimation  and  control,  none  of  the  large-scale  commercial  finite  element  systems 
has  facilities  for  error  estimation  or  adaptive  improvement. 

The  maturation  of  the  technology  of  estimation  and  control  of  discretization  errors, 
and  the  incorporation  of  this  technology  into  general-purpose  finite  element  systems,  will 
allow  the  analyst  to  select  only:  a)  the  initial  discrete  model  which  is  sufficient  to  resolve 
the  topology  of  the  structure  (or  the  geometry  of  the  domain);  and  b)  the  error  measure 
and  the  tolerance.  Then  the  finite  element  system  can  automatically  refine  the  model 
until  the  selected  error  measure  falls  below  the  prescribed  tolerance.  The  strategy  for 
adaptive  improvement  can  either  be  specified  by  the  user  or  automatically  selected  by  the 
program  (possibly  with  the  aid  of  an  Al-based  expert  system)  in  such  a  manner  as  to 


minimize  the  cost  of  the  analysis.  The  research  areas  that  are  needed  to  make  error 
estimation  and  adaptive  improvement  practical  are  identified  in  the  paper. 
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APPENDIX  I  -  ERROR  NORMS 

Error  norms  are  introduced  to  provide  quantitative  measures  of  the  magnitude  of 
the  error  in  the  finite  element  solutions.  These  are  measures  of  the  magnitude  of  the 
discretization  errors  in  certain  quantities,  and  therefore,  their  choice  depends  on  the 
goal  of  the  analysis.  Three  principal  choices  of  error  norms  are  commonly  used  in  finite 
element  methods:  the  energy  norm  ||  e  ||  ^  ,  the  mean-square  norm  !l  e  1  Iq  ,  and  the 
maximum  or  infinity  norm  ||  e  ||  <s>. 

The  energy  norm  is  defined  as  the  square  root  of  the  energy  of  the  error,  i.e. 

||  e  ||  =  (/  U (e)  dfi)1/2  ( 

E  n 

where  fi  represents  the  solution  domain;  U  represents  the  strain  energy;  e  is  the 
discretization  error  function  (difference  between  the  exact  and  the  finite  element 
solutions,  respectively). 

The  mean-square  (or  L  )  norm  measures  the  root-mean-square  of  the  errors  over 
the  solution  domain,  and  is  defined  by: 

||  e  ||  =  (/  e2  dfi)1/2  ( 
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The  maximum  norm  measures  the  maximum  absolute  value  of  the  error  over  its 


domain: 


e  =  max  e (x) 
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In  Eqs.  (A.2)  and  (A. 3),  e  refers  to  displacement,  stress  or  other  response  quantity. 

Of  the  three  error  norms,  the  energy  norm  is  the  most  commonly  used.  It  provides  a 
meaningful  quantitative  measure  for  the  error  in  the  finite  element  approximation  and  is 
equivalent  to  the  root-mean-square  of  error  in  stresses.  For  practical  engineering 
applications  the  maximum  norm  is  not  convenient  to  use. 


For  conforming  elements  based  on  minimization  principles  the  strain  energy  of  the 
error  decreases  monotonically  with  systematic  changes  of  the  discretization  (refining  the 
mesh  or  increasing  the  degree  of  the  polynomial  approximation). 

For  a  wide  variety  of  finite  elements  a  priori  error  estimates  have  been 
developed.  These  estimates  predict  the  asymptotic  rate  of  convergence  as  the  mesh  size 
tends  to  zero  (or  the  degree  of  polynomial  interpolation  function  goes  to  infinity). 

A  priori  estimates  for  stationary  (elliptic)  and  time-dependent  (parabolic  and 
hyperbolic)  problems  are  given  in  the  literature  on  mathematical  theory  of  finite 
elements.  Both  conforming  and  nonconforming  elements;  single-field  and  multifield 
models  have  been  studied. 

For  simplicity,  only  linear  elliptic,  self-adjoint,  positive-definite  problems 
(example,  linear  static  stress  analysis  problems  modeled  by  compatible  displacement 
models)  are  considered  herein.  For  these  problems  the  discretization  error  consists  of 
two  distinct  parts:  a)  approximation  error*,  and  b)  perturbation  error.  The  first 
(approximation)  error  characterizes  the  approximation  of  the  finite  element  solution  and 
depends  on: 

a)  Order  of  the  highest-order  derivative  used  in  the  formulation; 

b)  Degree  of  the  elements; 

c)  Dimension  of  the  largest  finite  element  (e.g.,  diameter  of  the  element);  and 

d)  Smoothness  of  the  solution. 

For  uniform  meshes  when  the  solution  is  not  smooth  (e.g.,  due  to  the  presence  of 
singularities  in  the  solution  inside  the  region  or  on  the  boundary)  the  convergence  rate 
decreases. 

Perturbation  errors  result  from  geometric  approximation  of  the  boundary; 
approximation  of  the  essential  boundary  and  continuity  conditions;  and  numerical 
approximation  of  integrals.  The  finite  element  solution  is  said  to  be  optimal  if  the 
perturbation  errors  do  not  exceed  the  approximation  errors. 


For  elliptic  problems  in  which  the  perturbation  error  is  small,  the  error  estimate 
can  be  given  in  the  following  simple  form: 

||  e  ||  <  Cf(N)  (A. 4) 

where  C  is  a  positive  constant  depending  upon  the  data  of  the  problem;  N  is  the  number 
of  degrees  of  freedom  in  the  model;  and  f(N)  is  a  function  of  N;  which  depends  also  on 
the  choice  of  interpolation  functions;  the  element  size;  and  the  smoothness  of  solutions. 

If  ||  e  ||  is  selected  to  be  the  energy  norm  II  e  |IE,  the  function  f(N)  can  be  given  by: 
f(N)  =  N_S  (A- 5) 

for  the  h  and  p  extensions  (Refs.  47  and  48);  and 

f(N)=l/eYN  (A. 6) 

for  the  h-p  extension  (Ref.  46). 

where  6  ,  y  and  0  are  positive  constants  which  can  be  determined  by  using  a  priori 
error  estimates. 

For  the  h-version  the  asymptotic  rate  of  convergence,  6  ,  is  the  absolute  value  of 
the  slope  of  the  plot  log  ||  e  ||  versus  log  N.  A  large  B  signifies  a  rapid  decrease  in  the 
error  with  increasing  N.  For  the  h-p  version  the  asymptotic  rate  of  convergence  can  be 
conveniently  obtained  by  plotting  log  ||  e  ||  versus  N9 . 


APPENDIX  II  -  HIERARCHIC  BASIS  FUNCTIONS 


Hierarchic  basis  functions  refer  to  nested  sets  of  shape  functions  which  have  the 
property  that  each  set  is  explicitly  contained  in  the  succeeding  sets.  The  nesting  of  the 
shape  functions  is  similar  to  that  of  the  terms  of  the  Fourier  series  expansion,  and  is 
usually  used  in  conjunction  with  the  p-extension.  However,  it  can  be  interpreted  in  the 
context  of  the  h-extension  when  the  shape  functions  associated  with  a  coarse  grid  are 
explicitly  contained  in  successively  refined  grids.  The  nesting  of  the  shape  functions 
results  in  the  nesting  of  the  elemental  stiffness  matrix  and  load  vector  (i.e.,  the  stiffness 
matrix  and  load  vector  associated  with  the  lower -degree  shape  functions  are  explicitly 
contained  in  successively  higher -order  approximations).  In  the  p-extension,  the  lower- 
degree  shape  functions  are  contained  in  successively  higher -degree  basis.  These  basis 
functions  are  constructed  from  integrals  of  Legendre  polynomials  and  are  referred  to  as 
p-hierarchic  shape  functions. 

For  one  dimensional  elements  and  two-dimensional  quadrilateral  master  elements 
with  C°  continuity  the  p-hierarchic  shape  functions  are  listed  subsequently  (see  Ref.  52). 
B.l  One-dimensional  Elements 

In  one  dimension  the  set  of  p-hierarchic  shape  functions  are: 

H1  =  1/2  (1-  O 
W2  =  1/2  (1+  O 

=  Jlizl  P>  (t)dt  ,  i  >  3 

l  2  j  1-2 

where  P-(t)  is  Legendre  polynomial  of  degree  i;  C  is  a  dimensionless  coordinate, 


The  shape  functions  have  the  following  orthogonality  property 


,  aw.  aw. 

/  ~aT  H1  dc  =  6ij 


i»j  >  3 


where  6  is  the  Kronecker  delta, 
ij 

Note  that  for  i  >  3,  W.  (-1)  =  W.  (+1)  =  0. 

i  i 

B.2  Two-Dimensional  Quadrilateral  Elements 


These  have  three  kinds  of  hierarchic  shape  functions: 

a)  Nodal  shape  functions.  These  are  the  4  standard  shape  functions  for  a 
quadrilateral  element: 

R  =  1/4  (l-O(l-n) 

W2  =  1/4  (l+c)(l-n) 

W3  =  1/4  (l+0(l  +  n) 

W4  =  1/4  (l-0(l  +  n) 

b)  Side  Shape  Functions.  There  are  4  (p-1)  shape  functions  associated  with  the 
sides  of  elements  (p  >  2).  The  shape  functions  associated  with  the  sides  n=-l  and  Z  =  +  l 
have  the  form: 


wj1}  =  1/2  (1 -n)  N.(?) 


1/2  (1+0  W.(n) 


(i =  3,  ...  p) 


with  similar  definitions  for  the  shape  functions  W^  and  W^  associated  with  the  sides 
n=  +  l  and  r=-l;  and  W  j(  z)  and  W  j(n);  i  >  3,  are  the  one-dimensional  shape  functions  defined 
in  the  preceding  subsection. 

c)  Internal  Modes.  There  are  1  / 2(p-2)(p- 3)  internal  modes  (p  >  4).  The  first  mode 

&[0)  (On)  =  (1  -Z2)(l  -  n2) 

The  other  internal  modes  are  obtained  by  multiplying  wj'0'1  by  Legendre 
polynomials  and  products  of  Legendre  polynomials,  as  follows: 


^0)  (C.n)  =  W^0)  P1(c) 

h\0)  U, n)  =  w{0)  P1(n) 

N<0)  (C.n)  =  _<0)  P2CC) 

N<0)  (c,n)  =  wj0)  PL(0  Px(n) 


with  similar  definitions  for  the  other  shape  functions. 

B.3  Two-Dimensional  Curved  Sided  Elements 

The  implementation  of  proper  mapping  procedures  is  important  in  order  to  ensure 
that  the  boundary  approximation  will  not  degrade  the  performance  of  the  element.  One 
approach  to  accomplish  the  mapping  is  through  the  use  of  linear  blending  function 
method  (see  Refs.  178  and  179),  in  which  the  Cartesian  coordinates  x,y  are  expressed  as 
functions  of  the  local  element  coordinates  ?  and  n  . 

The  inverse  mapping  can  be  obtained  by  using,  for  example,  Newton-Raphson 
technique.  Isoparametric  mapping  can  be  viewed  as  a  special  application  of  the  blending 
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a)  possibility  of  using  previous  solutions  and  computations;  and, 

b)  improved  conditioning  of  the  equation. 

3.  Hierarchical  shape  functions  have  been  developed  for  triangular  elements  with 
C1  continuity  (Ref.  ISO).  However,  the  form  of  these  functions  is  fairly  complicated  for 
practical  applications.  In  addition,  no  C*  shape  functions  exist  for  curved-sided 
elements. 
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